Stress Relaxation of Entangled Polymer Networks 

Gary S. Grest", Mathias Piitz^, Ralf Everaers'^, and Kurt Kremer'^ 
""Sandia National Laboratory, Albuquerque, New Mexico 87185 
Center for Microengineered Materials, University of New Mexico, Albuquerque, New Mexico 87106 
0^ \ '^Max Planck Institut fiir Polymerforschung, Postfach 3148, 55021 Mainz, Germany 

Q^ . The non-linear stress-strain relation for crosslinked polymer networks is studied using 

^~~' I molecular dynamics simulations. Previously we demonstrated the importance of trapped 

"t^ ■ entanglements in determining the elastic and relaxational properties of networks. Here 

^^ , we present new results for the stress versus strain for both dry and swollen networks. 

I Models which limit the fluctuations of the network strands like the tube model are shown 

,__! ' to describe the stress for both elongation and compression. For swollen networks, the 

. total modulus is found to decrease like (Vq/V)'^'^ and goes to the phantom model result 

only for short strand networks. 



O 



I. INTRODUCTION 



I I Understanding the elastic and relaxation properties of crosslinked networks and rubber has been a 

'^ longstanding problem [|l], |^ . One of the main unresolved issues has been the importance of trapped 

Q entanglements in determining the properties of the network. Edwards first introduced the idea that 

O obstacles produced by other chains give rise to a tube in which the monomers of a chain move. Later 

de Gennes generalized these same ideas to long chains, resulting in the reptation model for the dynamics 

^-H of an entangled polymer melt @,|§]. However how entanglements affect the properties of a crosslinked 

^ network remains controversial. The reason is twofold. The first is that most experiments do not allow 

\Q good control of all the microscopic parameters. Second, theoretical descriptions are very complicated due 

^^ to the presence of quenched disorder and usually contain several adjustable parameters, which are difficult 

2 to relate to microscopic details and are almost impossible to determine uniquely from experiment. 

Q>^ ' Computer simulations can directly address many fundamental questions regarding the dynamics of 

QN polymer networks. Details of the microscopic topology, such as the elastically active fraction of the network 

"H gind loop entanglements, can be identified and controlled. In particular, it is possible to isolate and quantify 

r^ their effects on macroscopic observables such as the elastic modulus. Simulation results on randomly 

I crosslinked, end-linked |^ and diamond networks [0-0| have conclusively demonstrated the importance of 

'^ trapped entanglements in determining the elastic and relaxational properties of the network. Results for 

Q the mean square displacement of the crosslinks show that they do not move completely freely in space as 

O assumed in the phantom network model but instead are confined to region of space of the order of the tube 

J> diameter of an uncrosslinked polymer melt 0. Further the monomers of strands of length A^^ > N^., where 

IT/ Ne is the entanglement length of the melt, are also confined to a tube. The fluctuations of the middle 

i-i monomers are found to diverge very slowly , as A^s in accordance with reptation theory. 

. - - , To test the theoretical models, it is advantageous to study model networks, preferably with equal strand 

lengths between crosslinks and no dangling ends. As such, we have carried out elongation simulations on 

ideal diamond lattice networks [0-^] and end-linked networks ^ . In the former, the network is constructed 

to have the coordination of a diamond lattice, with a well controlled topology. A second class of networks, 

which is not quite as ideal but experimentally realizable is end-linked networks [l^, 11|. In this case one 



starts with a melt of linear chains and adds /-functional crosslinkers to a fraction of the chain ends. Then 
by changing the temperature or other conditions, these groups can be activated, so that when another end 
comes in contact a chemical reaction can occur. We have made a number of model networks following this 
procedure [Q . A third class which is of high practical importance are randomly crosslinked networks, where 
crosslinks between melt chains are introduced between any pair of monomers. We have also made model 

1 



systems of this type, with the difference that we attached the end-monomers to a randomly chosen nearby 
monomer resulting in a tri-functional networks instead of the usual four-functional ones. The benefit of this 
method is that the networks contain no dangling ends which slow down relaxation and are not believed to 
contribute to zero-frequency stress. Here we present new results for the non-linear stress strain properties 
of end-linked networks as well as some results on swelling of randomly crosslinked networks. 

II. SIMULATION MODEL 

To simulate the networks we use the molecular dynamics method which has been successfully applied to 
study entanglement effects in polymer melts [|1^,|13[. In this model the polymers are represented as freely 
jointed bead-spring chains. All monomeric units of mass m interact via a purely repulsive Lennard- Jones 
potential 




m' + i 



r < Tc, 



Ujj{r) = { ''^[y-) -y-) ^4j ' ^'c. (1) 

r > Tc, 

where Vc = 2^/^cr. Monomer units connected along the chain or through the crosslinking procedure are 
connected by a finite extensible non-linear elastic (FENE) potential. The model parameter are the same 
as in ref. |]12[. The temperature T = e/ks- We use dimensionless units in which a = 1 and e = 1 and 



the basic unit of time r = cr(m/e)^'^. The simulations are carried out at a constant density p = 0.85cr~^. 
One of the most important features of this model is that the energy barrier for the crossing of two chain 
segments is high enough (r^ 70/cb^) such that crossing is virtually impossible. 

In our previous studies of the mean squared displacement of the chains in a dense melt, we found that 



the entanglement length N^, for this model is approximately 35 beads [|T2[. This result is obtained by 
comparing the crossover time t^ at which the motion of the monomers in the chain slows down from its 
initial Rouse-like motion at early times to the slower t^^* behavior predicted by reptation theory. This 



result was confirmed by come recent simulations [|T3[ on very long chains (A^ = 350,700 and 10000) which 



also found that the segment motion of the chains was in very good agreement with the reptation model. 
Estimates from these simulations suggest Ne = 32 ± 2 [|^. However calculations of the plateau modulus 
G%, which were not carried out in our earlier studies due to lack of computer time, indicate that Ne is 
much larger. This estimate of Ne^p was obtained using the standard formula of Doi 

_ A pkeT 

^" - 5 A^- ^'^ 

Depending on the formula we use to fit the non-linear stress strain data, we obtain values for N^^p between 
65 to 80, about twice as large as that obtained from the displacement of the chains. We use the average 
value 72 when comparing our results to experiment. This suggests an error in the pre-factor | in Eq. |[ 
Experimentally A'e is only determined from G°^. We also find that diffusion constant data from the 
simulations and from experiment agree very well when both are scaled by their Rouse value for short 
chains and plotted versus N/N^^p |T^ 



Here we present results for end-linked systems (/ = 4) consisting of M chains of length A^, where 
M/N = 1000/35, 500/100, and 120/350. Due to the slow relaxation and long simulation runs required to 
obtain the stress at a given strain, the results presented here are for one configuration for A^ = 35 and 100 
and two configurations for A^ = 350. For comparison we also present results for the plateau modulus for 



A^ = 350 and A^ = 700 |T^. For details of the end crosslinking procedure, see ref. [0. For the swelling 
studies, we present results for randomly crosslinked systems (/ = 3) of size M/N/Ng = 500/700/233 and 
1600/50/18, where A^ is the melt precursor chain length before crosslinking and A^^ = A^/3 is the average 
length of a network strands between crosslinks. We will restrict our analysis to the effect of (partial) 
swelling in good solvent on the stress-strain relationship of these networks. 



III. NON-LINEAR STRESS STRAIN IN THE DRY STATE 

The elastic modulus and non-linear stress versus strain are easily measured experimentally for networks. 
Numerically the shear modulus can be obtained from uniaxial, volume conserving elongation of the sample 
which can be described by a diagonal deformation tensor A with the stretching factor A^^; = A and the 
contraction factors Xyy = Azz = A~^/^ in the other two directions. The normal stress a^ is then readily 
determined from the microscopic virial tensor, a^ = (Jxx — ipyy + cr2^)/2, where x is the direction of 
elongation. In the present simulations we vary A from 2.5 (extension) to 0.5 (compression). Due to finite 
system size, we are limited to relatively large strains A > 1.2 or A < 0.8. For smaller strains, the stress is 
small and difficult to determine from the noise. After an initial step strain, the strain decays very slowly. 
For the longest chain lengths studied (A^^ = 350), runs typically of at least lO^r (< 2 x lO^r for shorter 
network strands) were needed to reach equilibrium. 

The classical theory of rubber elasticity |jl| is based on the empirical fact that the network strands have 
Gaussian coil statistics. With the exterior deformation described by the tensor A the resulting free energy 
stored in the network strands (of entirely entropic origin) is: 

where V is the volume of the sample, Pchain the number density of elastically active strands, / their 
functionality and Ka are the diagonal elements of the deformation tensor. The above result was derived by 



James and Guth [Tj] under the assumptions that the crosslinks and the connecting strands are allowed to 
fluctuate freely without mutual hindrance of each other (this is commonly referred to as the phantom chain 
model). If one fixes the crosslinks in space so that they deform affinely under strain (the affine model) the 



factor (1 — 2//) disappears and on arrives at the old result given by Kuhn [T^. With this on can readily 



derive an expression for the normal tension in the sample as function of the elongation [16 



— ^^ - "-^^^^ - 7) (^^ 4 1 • w 



where L^/A is the transversal cross section of the sample. Flory and Erman [17, 18| extended the theory to 
allow entanglement effects on the crosslinks by reducing their fluctuation radii, thus interpolating between 
the affine and the phantom model with a more complicated stress strain behavior. The model is quite 
successful describing many experimental stress-strain relations as well as our data (see Fig. |l]a). However, 
due to our additional knowledge on the conformations of the chain ||^ , we know that the fluctuations of 
the entire chains are restricted and therefore we think that constraint models for junctions alone should be 
discarded. Many models exist which limit the fluctuations of the network strands like the tube-model by 
Heinrich et al. [|l^ or some recent constraint models by Kloczkowski et al. ||20|, Rubinstein and Panyukov 
(RP) [pli[2^ and Everaers [^. Of these, especially the RP model is rather appealing, since it introduces 



only one additional parameter, the chain length between entanglements N^- Erman's and Everaers' recent 
models introduce a large set of constraint parameters and little is known yet how to choose them, although 
the physical meaning of these parameters is relatively clear. At least for Everaers' model it has been 
demonstrated how to obtain more detailed knowledge about network fluctuations by measuring these 
parameters directly through simulations of highly idealized networks with diamond lattice topology P] . 
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FIG. 1. Stress-strain relation for end crosslinked networks compared to (a) the Mooney-Rivlin, Eq. [^ (solid 
lines), and constrained junction (non-affine deformation parameter k = 5.0) fl^ (dashed lines) models and (b) 
the Rubinstein-Panyukov form, Eq. (g). Data are for A'^, = 35 (□), 100 (o) and 350 (triangles). Also shown is the 
plateau value for an uncrosslinked melt of chains of length A^ = 700 (■) . 



An analysis of the type presented in [^] is beyond the scope of this paper and we shall limit our analysis 
along the lines of the RP model which has an extra entanglement term in the free energy: 



i=l ^ ' 



(5) 



It is customary to renormalize the stress a^ by the classical functional dependence on A, a*^ = o"Ar/(A^ — 1/A) 



to emphasize deviations from this behavior. Together with some recent improvement p2[ on the theory 
which takes into account reorganization of the constraining tubes due to anisotropic tube deformations |^ 
this additional term results in a stress-strain relation for uniaxial strain: 



a^ = G, + 1.84Ge/(A + OMX~^/^] 



(6) 



In Fig. ^ we show the stress-strain curves for our various model end-linked networks together with fits to 
Eq. P to the semi-empirical Mooney-Rivlin [|l| formula 



a^ = 2Ci + 2C2/A 



(7) 



and to the resulting stress-strain relation of the theory of Flory and Erman (formula not shown due to 
complexity) for comparison. We can see that the elongation branch can be described by each of these 
equations, while the Mooney-Rivlin form provides a clearly inadequate description of the compression 
data. Note, that the stress has been normalized by the classical stress dependence (Eq. ^ on A, so that 
deviations from the classical behavior are more clearly visible. 

However a more quantitative analysis of the resulting fit-parameters shows that the situation is less 
conclusive. We performed two fits to the PR model. For the first fit the value of Ge was determined from 
the plateau modulus for the uncrosslinked melt of free chains and was held fixed. Using the RP form, 
Ge = 0.0102e/(T-3, while G^ = 0.0145, 0.0097, and 0.0034e/or"3 for iV, = 35, 100 and 350, respectively. If 



Gc is equavalent to either phantom or affine results the values of Gc{Ns) should scale with the crosslink 
density (N-^). The values for the ratios Ge(35)/Gc(350) = 4.3 and G'c(100)/Gc(350) = 2.9, compared 
to the phantom/affine prediction of 10 and 3.5. The absolute values for Gc generally appear too high 
compared to the phantom and/or the affine results. A different way to look at the results is to set Gc to 
their expected values and fit Gc- For our shortest strand network A^^ = 35, fluctuations of the crosslinks 
are of the order the radius of gyration of the entire chain, hence we expect Gc should equal its value for the 
phantom model, Gc = 0.012e/a~^. Using this value for Gc, we find Ge = 0.0014e/cr~^. For our Ns = 350 
networks the crosslinks explore much less space than the radius of gyration of the whole chain, thus the 
behavior is expected to be more in accordance with the affine model, Gc = 0.0022e/cr~^ (after adjusting for 
the fact that 10 of 120 chains are not elastically active). Using this value for Gc we find Ge = 0.0012e/cr~^. 
Thus the effective Ge becomes larger as the chains become shorter, i.e. the presence of more crosslinks 
provides stronger constraints as in the melt case. A more rigorous test of the RP model would be to 
compare networks of strand lengths Ns = 350 and 700, which is beyond the amount of available computer 
time at present. However experimentally the more relevant region in parameter space lies in the range 
A^ = 35 to A^ = 350 since it is difficult to make networks with very large strand lengths. 
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FIG. 2. Plot of modulus for end linked networks from simulation (o) and experiment (A). The modulus has 
been normalized by its value at Ne- The experimental data for G°^{Ms) of PDMS networks are from Patel et 



ai. in|]. 



The elastic modulus G°^ can be obtained from extrapolating the slope of the non-linear stress-strain 
curves to A = 1. Results for G"^ are shown in Fig. |l| using the RP form to extrapolate to A = 1. The 
MR form gives larger values for G°^. To compare to experimental data for PDMS end-linked networks, we 
have normalized the modulus for its value at Nc^p = 72. For PDMS the entanglement molecular weight 
Me,p = 9600. As seen from Fig. |l|, our simulation results agree with experiment. The large scatter in the 
experimental data for strand lengths Ng '^ 2 ■ Ne^p shows the difficulty to create good networks with the 
endlinking method. 



IV. SWOLLEN NETWORKS 



Experimental results [P1,P3| show that the entanglement contribution to the modulus decreases signifi- 
cantly if the networks are swollen in good solvent. To investigate this effect we let our randomly crosslinked 
networks swell in vacuum. Since our model only involves repulsive interactions vacuum acts as a perfect 
solvent for our model networks. At various degrees of (incomplete) swelling we performed stress-strain 
experiments (only elongation) with two of our swollen samples. We used formula (^) to fit the resulting set 
of stress-strain curves and estimate the modulus, though for extracting the modulus any of the forms for a^ 
could have been used. Fig. ^ shows the stress-strain relations for the unswoUen states and at equilibrium 
swelling (when elastic free energy and entropy of dilution cancel each other). The swollen systems show 
the empirically expected behavior of classical rubber elasticity (^ which becomes clearer when one looks 
at the classical and entanglement contributions to the modulus shown in Tab. |. These were determined 
by fits to Eq. (P). One can see from the table as well as from Fig ^. that the total modulus (slope of the 
curves for A ^ 1) approximately decreases like g~^ = {Vq/VY'^. For the short strand network Ng = 18 it 
is nice to see that the system recovers the phantom result at maximum swelling (the effective values for Gc 
and Ge are no longer meaningful, since finite chain length effects become very important at high swelling 
and large deformations which are not accounted for in the theory, but still they are useful to extract a 
total modulus). For the long chain a entanglement contribution remains, however it is shifted over into 
Gc which is still about a factor of 1.5 to 4 larger than the afiine and phantom values, respectively). If 
RP theory is applied to swelling, in the spirit of Flory and co-workers (i.e. assuming that the statistics 
of the chains remains Gaussian with swelling), this theory (as well as others) actually predicts, that the 
entanglement contribution should vanish and the classical crosslink contribution should remain constant 
which qualitatively agrees with our findings. 
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FIG. 3. Stress-strain relations for a randomly crosslinked three-functional system of average strand length 
Ns = 233 for various degrees of swelling: q = 1.0 (•), q = 1.2 (0), q = 1.4 (a), 1 = 1.6 (□), q = 1.9 (o). For clarity 
the curves are shifted in y-direction by O.Ole/cj^ with respect to each other. 



This follows from a simple analysis of the leading powers in A of the different contributions to the free 
energy. Let us split up the deformation into an isotropic swelling factor q and an anisotropic deformation 
A^j : Ajj = qh.'ii- The free energy due to crosslinks is then proportional to g^, but the system size L also 
changes proportional to q. Therefore the powers of q in the formula for the stress cancel exactly. The 
entanglement contribution in the RP model is proportional to the first power of q only and thus its stress 
contribution should decrease like 1/g, whereas we find it behaves more like g~^ for our networks. The 
most important key to understand this behavior certainly lies in the non-Gaussian statistics of the swollen 
chains. We recently found, that the fractal structure of the chains changes significantly under swelling and 
can be described as random walks only on length scales beyond the entanglement length Ne and below can 
be described by a fractal exponent of u 



0.7 26 



which is in disagreement with the widely assumed single 
chain good solvent exponent v ~ 3/5. Further our results for long-chain networks seem to favor a form of 
the free energy where the entanglement contribution behaves approximately like Fe{q., A) -^ G{Ns) ■ Fc{X) 
for large q. Loosely speaking, the trapped entanglements behave more like physcial crosslinks for high 
swelling. 



Ns 


Gph 


Q 


Gc 


Ge 


Gtot 






1.0 


0.0102(05) 


0.0110(10) 


0.0212(15) 






1.2 


0.0074(06) 


0.0058(10) 


0.0132(16) 


233 


0.0012 


1.4 


0.0062(05) 


0.0037(02) 


0.0194(12) 






1.6 


0.0059(03) 


0.0018(04) 


0.0077(07) 






1.9 


0.0058(05) 


-0.0005(05) 


0.0053(10) 






1.0 


0.0190(20) 


0.0120(20) 


0.0310(40) 






1.2 


0.0254(60) 


-0.0047(60) 


0.021(12) 


18 


0.017 


1.3 


0.0327(30) 


-0.0166(30) 


0.0161(60) 






1.4 


0.0314(12) 


-0.0162(15) 


0.0151(27) 



TABLE L Classical crosslink (Gc) and entanglement (Ge) contributions to the total modulus Gtot in swollen 
(g^ = V/Vo) randomly crosslinked three-functional systems of average strand length Ng. Here Gph = Pchain{^ — '^/f) 
with / = 3 is the modulus predicted by the phantom model (y). Errors in parenthesis are those obtained by the 
non- linear fit. 



V. CONCLUSION 

In this paper we have presented molecular dynamics simulation results for the non-linear stress strain 
properties of a crosslinked polymer network. Because the topology of the network can be readily con- 
trolled in the simulations, computer simulations provide a powerful way to test various theoretical models. 
Combined with our earlier studies, we are now beginning to obtain a better microscopic understanding 
of how the topology of the network controls the dynamics and relaxation of the chain segments and the 
crosslinks. All of our data strongly support constraint models in which the monomers are constrained to 
move in a tube formed by the neighboring chains. The crosslinking process traps in the dynamic entan- 
glements present in the uncrosslinked melt. These trapped entanglements dominate the stress relaxation 
and modulus when the strand length between crosslinks is large. In the unswollen, dry state, the chains 
are Gaussian in agreement with all models. However as the system is swollen in a good solvent, they 
become strongly non-Gaussian on length scales less than the distance between crosslinks. In the dry state, 
a number of model for the stress-strain relation can be used to describe the data, although further work 
is necessary to arrive at a more quantitative understanding physical of the resulting fitting parameters. 
Clearly a more microscopic approach investigating network fluctuations is necessary to distinguish between 
these models. In the swollen state we find that short stranded networks eventually recover the phantom 
result for the modulus, however finite chain length corrections and a proper account for the non-Gaussian 
statistics are clearly needed to descibe their entire (A 7^ 1) stress-strain relations. In long strand networks 
the entanglement contribution to the modulus decreases with swelling but compared to the predictions of 
the phantom and affine models a significant contribution remains. The stress strain curves at high swelling 
are better described by classical (crosslinks only) theories, hence the remaining effect of entanglements 
appears more like an additional crosslink contribution. 
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